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Most studies of GeV Galactic Cosmic Rays (GCR) nuclei assume a steady 



Ph| state/continuous distribution for the sources of cosmic rays, but this distribution 

O ' is actually discrete in time and in space. The current progress in our understand- 

^ [ ing of cosmic ray physics (acceleration, propagation), the required consistency in 

^ ' explaining several GCRs manifestation (nuclei, 7,...) as well as the precision of 

present and future space missions (e.g. INTEGRAL, AMS, AGILE, GLAST) 
^ ' point towards the necessity to go beyond this approximation. A steady state 

I semi-analytical model that describes well many nuclei data has been developed 

in the past years based on this approximation, as well as others. We wish to 
extend it to a time dependent version, including discrete sources. As a first step, 
the validity of several approximations of the model we use are checked to validate 
the approach: i) the effect of the radial variation of the interstellar gas density is 
inspected and ii) the effect of a specific modeling for the galactic wind (linear vs 
constant) is discussed. In a second step, the approximation of using continuous 
sources in space is considered. This is completed by a study of time discreteness 
through the time-dependent version of the propagation equation. A new analyt- 
ical solution of this equation for instantaneous point-like sources, including the 
effect of escape, galactic wind and spallation, is presented. Application of time 
and space discretness to definite propagation conditions and realistic distribu- 
tions of sources will be presented in a future paper. 
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Subject headings: Cosmic rays 



1. 



Introduction 



The cosmic ray flux at any given position in the Galaxy is due to many sources, which 
are probably related to the remnants of supernovae. Each source of position rj and age 
yields, at position tq and time t = (now), a flux N{fi,ti,fQ,tQ = 0) that can be obtained 
by solving the diffusion equation with the appropriate source term and boundary conditions 
(see below). The total flux is given by 



This model has been coined Myriad model by Higdon and Lingenfelter (2003). 

The effects of discreteness have been studied in the past. As regards the spatial dis- 
creteness, Lezniak and Webber (1979) start with the time dependent diffusion equation when 
both spallations and energy losses are taken into account, to eventually derive the steady- 
state Green function necessary to study the no-near source effect (expected to reproduce a 
depletion of the path length distribution at low grammage). As for the temporal discreteness, 
Owens (1976) derived time dependent solutions in the framework of halo models (see also 
Lezniak 1979), including spallations but assuming a gas density which is constant throughout 
the diffusive volume (in this case the mean time and the mean matter crossed are propor- 
tional). The most complete work to provide "simple" formula to the time- dependent case is 
probably Freedman et al. (1980), who derive the mean age and the grammage distribution 
to seek if it is possible to constrain the propagation parameters from current observations of 
charged nuclei. They conclude that a large degeneracy in propagation parameters remains 
in most cases. It has been shown that the distribution of electrons and positrons at high 
energy is particularly sensitive to nearby sources, due to their huge energy losses (Aharonian 
et al. 1995, see also Coswik and Lee 1979). A nearby source such as the supernova remnant 
(SNR) associated with the Geminga pulsars probably has a great influence on them. The 
situation is less clear about the stable charged cosmic ray spectra, but nearby sources are 
expected to be more and more important as energy grows. Whereas some authors estimate 
the contribution from Geminga to be at most 10% (Johnson 1994), some others (Erlykin and 
Wolfendale 2001a) argue that a single supernova event can explain the feature observed in 
the spectrum at a few PeV. This could be the Sco-Cen association that is expected to be able 
to generate the Local Bubble 11 Myr ago (Benftez et al. 2002). This can be tested by mea- 
suring the anisotropy (see e.g. Dorman et al. 1985), or by gathering observations from the 
past and from the rest of the Galaxy. For instance, Ramadurai (1993) argued that Geminga 
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could be responsible for an increase of the Cosmic ray flux by a factor 1.8 by inspection of 
Antarctic ice sediments. The measurement of the isotopic composition of the Earth crust 
(see e.g. Knie et al. 1999), of meteorites, and of ice cores may be used to investigate the 
time variation of the cosmic ray flux on Earth (see also Erlykin and Wolfendale 2001b). 

Most studies of the chemical composition of cosmic rays assume that the sum (1) may 
be approximated by an integral 



which is equivalent to (1) only in the limit of a source distribution which is continuous in 
space and time (continuous distributions for rj and tj). This approximation is justified if the 
sources are numerous and densely distributed, but it probably fails for nearby and/or recent 
sources, for which the detailed location and age should be known. This paper is devoted to 
investigate the validity of this approximation and to provide a more accurate description of 
diffusion when it fails. The chemical composition of cosmic rays is determined by the quantity 
of matter that has been crossed by primaries during their propagation from the sources to 
the Earth, and it can be conveniently described and studied by the grammage distribution. 
Sec. 2 is devoted to the study of the grammage distribution in diffusion models, and the 
importance of space discreteness of the sources on the cosmic ray composition is investigated. 
Sec. 3 presents a new analytical solution of the time-dependent diffusion equation for point 
sources, taking into account spallations, convective wind and escape. It is then used to find 
a criterion to separate the sources in two categories, one containing the faraway and old, 
which can be modelled by the usual steady-state model, the other containing the close or 
recent, which require a finer description. Sec. 4 concludes and presents some applications of 
the present work which will be further developped in a next paper. 



During their journey between a source and Earth, Cosmic Rays cross regions in which 
interstellar matter is present. The nuclear reactions (spallations) induced by the collisions 
lead to a change in the chemical composition. Cosmic Rays can be sorted according to the 
grammage, i.e. the column density of matter they have crossed (denoted by x in this paper 
and usually expressed in g cm~^). If it is temporarily assumed that Cosmic Rays do not inter- 
act with the matter they cross, i.e. the spallations are switched off, their density originating 
from a source located at fg, detected at r*o, and having crossed the grammage x is called 
the path-length distribution. It can be used to compute the probability of nuclear reaction 
when the spallations are switched on, and thus it provides a tool to compare several diffusion 
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models, as similar path-length distributions give rise to similar chemical compositions. 



2.1. Definition - Generalized diff'usion equation 

The evolution of the path-length distribution Q{f,t,x) at position r and time t can be 
described by a generalized diffusion equation inspired by Jones (1979): 

^^i^lA = KAGir, t, x) + q{r, t)S{x) - ^^^^ll^rnvnU^I , (3) 

where in the right-hand side, the first term stands for diffusion, the second is the source 
term (creating particles with null grammage) and the third gives the augmentation of the 
grammage due to the crossing of matter. It should be noted that when the density nism is not 
homogeneous in the whole diffusive volume, grammage is not proportional to time (this will 
be discussed further in Sec. 2.2.3). We emphasize that Eq. (3) does not take into account 
the influence of the spallations on propagation, but rather introduces the variable a; as a 
counter which keep tracks of the quantity of matter crossed by the Cosmic Rays, spallations 
being switched off. 

The density can finally be obtained from the grammage distribution by noting that 
a primary cosmic ray having crossed a grammage x has a survival probability given by 
exp(— (Tx/m), where a is the destruction cross-section. In the following, m will denote the 
mean mass of the interstellar medium atoms. The probability for a Cosmic Ray emitted in 
fg and reaching fg unharmed is written as 

poo r (T ~i 

Np{fs, To) = expl X \ g{fs, fo, x) dx . (4) 

Jo I m J 

For a secondary species, a similar expression can be written 

POO 

Ns{fs,fo) = / gs{x)Q{fs,fo,x) dx (5) 







where the function gg is obtained by solving the set of equations 

-1- = 9s + > Qi (6) 

dx m ^-^ m 

with the initial conditions (the values of gs for x = 0) set to the source abundance of the 
considered species. 

The path- length distribution, along with the B/C ratio, are computed in the next para- 
graph. We take the opportunity of this computation to consider a more general situation 
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than in our previous works, namely by considering a realistic radial dependence of the matter 
density in the disk. The aim is two folds : first, we want to comment on the computation 
of path lengths by Higdon and Lingenfelter (2003), and in particular we want to discuss 
the existence of a feature at x ~ 2 g/cm^ which they claim is present because of the H2 
ring in the galactic disk. Second, this provides a way to effectively take into account the 
radial dependence by computing the mean matter density which is probed by each cosmic 
ray species, and to introduce this mean density back in our code. 



2.2. Analytical result for radial distribution of matter 

2.2.1. Path Length Distribution (PLD) 

We compute the grammage distribution in the diffusion model we used in previous 
studies (Maurin et al. 2001; Donato et al. 2002). It exhibits cylindrical symmetry, escape 
happens through the z = ±L and R = 20 kpc boundaries, galactic wind is constant in the 
halo and matter is localised in a thin disk at 2; = 0. We consider a radial dependence of 
the surface mass density S(r) taking into account the radial distributions of HI, HII and 
H2 (Ferriere 1998). For convenience, we normalize this quantity by the local mean surface 
density = 6.2 x 10^° cm~^, i.e. we introduce /(r) = S(r)/S?j^. The models used in our 
previous studies considered only flat matter distribution, to keep the problem tractable in a 
semi-analytical way, and we take the opportunity of this study to investigate the importance 
of this assumption. 

The generalized diffusion equation (3), with the left-hand side set to 0, can then be 
solved as detailed in Appendix B, by expanding the quantities over a set of Bessel functions. 
The solution is 

^(-■-■-'-f: ^""^^^T' -M^l)E-..^---eW, (7) 

1=0 ' ^ j=0 

where 9(x) is the Heaviside distribution, the Uij and the Xj are the eigenvectors and eigen- 
values of the matrix 

= KS-Jf{Q) ^^^^ (^) MQp)MQp)f{pr) dp with ^^ = ^, 
For a fiat distribution of matter (/ independant of r), this expression reduces to 

mwE? , / SiL\ ^ 
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We have not considered energy losses in this computation, as the aim is not to provide a 
very sophisticated modeUing of the cosmic rays diffusion, but rather to give an estimate of 
various effects. It follows that the results presented here do not apply directly to electrons 
and positrons, for which the energy losses are predominant. 

2.2.2. Application 

We now use the above expression to evaluate the effect of the choice of S(r) on the 
path-length distribution, and then on the composition of cosmic rays. This is done by first 
computing the PLD for a flat and for a more realistic matter distribution. The result for 
K = 0.03 kpc^/Myr and L = 5 kpc is displayed in figure (1). The feature at x ~ 2 g cm~^, 
visible in the Higdon points (crosses) is never reproduced by the analytical result. This 
difference is discussed in the next paragraph. 

The effect on the composition of cosmic rays is illustrated in Fig. 2, where the B/C ratio 
is computed in the two cases from the PLD. For low values of the diffusion coefficient Ko, 
this ratio is not very sensitive to the global distribution of matter, as the diffusion range is 
smaller, wheras for higher values of Kq, the B/C ratio is sensitive to the increase of matter 
density at r ~ 4 kpc. As the diffusion coefficient K{E) actually increases with energy, the 
spectrum is likely to be affected by this effect, cosmic rays of higher energy probing a larger 
portion of the galactic disk. However, for high values of i^, the difference between a fiat 
and a realistic distribution becomes independent of K, hence of i?, so that the shape of the 
spectrum is not affected. In particular, if Kq is high, as hinted at by Ptuskin and Soutoul 
(1998), then the spectrum is not affected by the presence of a H2 ring. 

These results provide an effective way to take into account the radial distribution of 
matter in our model. For a given energy, we can compute the constant interstellar gas 
density n"!^^^^ that must be assumed to reproduce the observed B/C ratio. The same thing 
can be done for sub-Fe/Fe, providing another effective density ri'H^\ These densities are 
different because the corresponding species have different diffusion ranges, the latter being 
much more sensitive to spallations. The flux of each species can then be computed using the 
appropriate effective interstellar matter density. 

2.2.3. Comparison to the Higdon and Lingenfelter approach 

Another approach has recently been proposed in Higdon and Lingenfelter (2003) to 
compute the grammage distribution: starting from the CR distribution w{f,t) due to an 
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instantaneous point source of age t, the mean matter density seen at age t by these CR is 
computed as 

// / d^fw{f,t) nisrr,{r) 

Usnit) = — 




d^rw{r, t) 

Summing all the individual contribution with the appropriate weight, the mean matter 
density ricrit) seen by the CR of age t reaching the Earth is computed. From this quantity, 
the evolution of the mean grammage is computed as a function of time, through 

dX 

— = rhpcucrit) (9) 

which corresponds to their Eq. (4.1) and (4.2). We argue that this approach is not correct, 
for the following reasons. 

First, the derivation of the time evolution of X{t) is given in Appendix B, and the 
equation (9) is not recovered. The above expression would be correct only with another 
(tricky) definition of ncr{t) and nsn{t)- 

Second, the averaging process (8) gives the mean density as seen by all the CRs emitted 
by the source, whereas what is needed would be the mean density seen by the CRs that 
reach the Earth (those we do observe). These quantities are different, and the corresponding 
time evolution is computed in Appendix B. 

Finally, and most importantly, it is quite tricky to infer the grammage distribution 
Q{f,x) (needed to apply the weighted slab technique) from the mean grammage X(t), and 
their relation (4.3) is not correct. To see that more clearly, consider the more fundamental 
quantity Q{f,x,t) giving, at position r, the density of CRs having crossed a grammage x 
at age t. The quantities introduced by Higdon and Lingenfelter (2003) are then related to 
Q{f, X, t) through 



dxxQ (r , X, t) 

w{f,t)= I dxg{r,x,t), g{f,x)= I dtg{f,x,t), x{t) = ^^ (lo) 

dx Q (r , x, t) 

which are fundamentally different from their (4.2) and (4.3). In particular, their X{t) is ac- 
tually X{t) and should be a function of t in their (4.1) and (4.2), whereas the X that appears 
in their grammage distribution (4.3) should be a parameter and as such is independent on t. 

This explains why we do not find the same grammage distributions as Higdon and 
Lingenfelter (2003). Their assimilation of the grammage distribution from individual sources 
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to Dirac distributions has the effect of sharpening the final grammage distribution. In 
particular, the feature at x ~ 2 g/cw? is not present in our results. The distributions appear 
to be actually quite close to exponentials., i.e. to Leaky Box distributions. 



2.3. Spatial discreteness of the sources in a steady-state model 

2.3.1. General results 

We now want to investigate the effect of discreteness of the source distribution on 
the Cosmic Ray composition, through the path-length distribution. We first compute this 
quantity for a point source, and we then compare the path-length distributions obtained 
for a set of point sources and an equivalent continuous source distribution. For the sake of 
simplicity, we focus on the case of a uniform distribution of matter, for which Aij is diagonal 
and the solution given in App. B can be simplified. As the composition of Cosmic Rays is 
only measured in the galactic disk, we express the results in 2; = 0. It is then found that 

Qi{z = 0,x) = Q exp {-x/xi) Q{x) 

with 

x. = ^^('^ + ^coth^^^^ 



2K V^w 2 V 2 

and r„ = 2K/Vc. This expression could have been obtained by an inverse Laplace transform 
of Fourier-Bessel coefficients of the steady-state density (see e.g. Maurin et al. 2001) 

iV,(^ = 0) '^^ 



vY^^^a + K + KSi coth {SiL/2) 



The Qi are obtained by Fourier-Bessel transforming the radial source distribution, which 
is assumed to be point-like and located in the galactic disk. Unless this point source is at the 
galactic center, the cylindrical symmetry of the problem is broken and the previous study 
does not apply. However, as the influence of the R boundary is expected to be negligible, 
we consider that the diffusion volume is not limited in the radial direction [R 00). The 
origin can then be set at the position of the source, which restores cylindrical symmetry. In 
this R ^ 00 limit, the summations over Bessel functions become integrals, the discrete sets 
gi, Xi, Ni become functions q{k), x{k), N{k) and the final result is obtained by performing 
the substitution 1/Jf{(i) — + k7rR/2, J2i ^ I d^Rk/n) and Q/R k, 

^(r,2; = 0,x) = — 4— / dk Jo{kr) q{k) exp [ ^) (11) 
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with 



l/r„ + S{k) coth{S{k)L) 
and 

= 2 (1 + A:^y'' , q{k) = kj\dr^-^^Mkr) = ^ 
In the particular case of infinite L and r„ = 0, the expression (11) gives 

1 K^x ( x^K^ 



which is the expression obtained in Taillet & Maurin (2003) from a random walk approach. 



2.3.2. Impact on the chemical composition 

The path- length distributions (11) are displayed in Figs. 3 and 4, for several propagation 
conditions, to emphasize the relative effect of escape, galactic wind and spallation. The 
effect of escape or wind is more important at higher grammages. This was expected as 
they correspond to longer paths. The effect of the convective wind is seen to be similar to 
that of escape, but quantitatively different at low grammages. To understand this, let us 
first consider diffusion in free space without wind. Several kinds of paths are responsible 
for low grammages: short paths, connecting us to nearby sources in the disk, and longer 
paths that wander in the halo without crossing the disk too much. It happens that the 
second kind is rather important, which explains why escape from a boundary at 2; = ±L 
(which kills the paths wandering too far in the halo) actually affects the low end of the 
grammage distribution. When wind is present, the short paths are more important, and the 
low grammages are less affected. 

Even for a point source, the path-lengths are quite broadly distributed around the mean 
value, as can be seen in Figs. 3 and 4. As a result, the grammage distribution due to a set of 
discrete sources is smoothed to a great extent, and that it is very unlikely to have observable 
consequences on the composition of Cosmic Rays on Earth, except for very nearby sources 
(see below). To illustrate this point, we compare the grammage distributions and the B/C 
ratio from a smooth distribution to that obtained from a discrete sample representative 
of this distribution. The relative importance of the sources located at different distances 
in making the observed composition (e.g. the B/C ratio) was estimated and discussed in 
Taillet and Maurin (2003), and it was found that the nearby sources can be responsible for 
a substantial fraction of the flux of each species. The contribution of a point source located 
at a distance r to the B and C flux, obtained from (11), is displayed as a function of r 
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in Fig. 11. More specifically, we divide the sources into the far (r > L), the intermediate 
(O.IL < r < L) and the nearby (r < O.IL), and Table 1 gives the fraction of flux coming 
from these regions, for different species. The grammage distribution from a ring dehmited 





B C sub-Fe Fe 


'^B/clkpc 2) 


< r/L < 0.1 
0.1 < r/L < 1 

1 < r/L 


3 % 12 % 7 % 20 % 
57 % 65 % 71 % 68 % 
40 % 23 % 22 % 11 % 


100 
20 



Table 1: The first four columns give, for the four species indicated (B, C, sub-Fe and Fe), 
the fraction of flux due to the three rings defined in the left column (0 < r/L < 0.1, 
0.1 < r/L < 1 and 1 < r/L), in the case of a homogeneous source distribution. The last 
column gives ub/c, the surface number density of discrete sources (in kpc~^), randomly 
distributed in these rings, beyond which the difference in flux with the continuous case is 
smaller than 1 %. There is no number given for the inner ring because there is no motivation 
to randomly draw the sources located very close to us, as these should be observed. The 
quantity ub/c was computed for L = 1 kpc and should be divided by L^ for other values. It 
does not depend on the value of K. The effect of the convective wind can be estimated by 
replacing L hj L* = K/Vc = r^/2. 



by i?min and i?max is given by 

2^SP Rl^^ - R 



g^r,z = 0,x) = ,_^.to ^ 02 (12) 



ism "max "min 

00 



dk {R^a^Ji{kR 

max ) Rm'm 

We randomly draw the position of sources inside each ring, with a given surface density of 
sources, we compute the corresponding flux, and the relative difference with the smooth case 
is computed. As expected, this difference gets smaller as the source density is increased, 
i.e. as the granularity of the source distribution is decreased. Table 1 gives ub/c, the 
surface number density of sources in the disk (in kpc~^) beyond which the difference in 
composition between the discrete and continuous case is less than 1 %. The secondaries 
are less sensitive to the discreteness of sources, as their sources (the primaries) do have a 
continuous distribution. These results do not depend on the value of K, as it is only sensitive 
to the relative importance of nearby and remote sources. 

The results are presented as a function of r/L. If we now consider the disk defined 
by r < 1 kpc, the resulting effect of granularity is quite insensitive to the value of L (or 
Vc), as long as L > 1 kpc (or r^ > 1 kpc). This is because for r < L (or r < r^,), the 
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cut-off effect of escape (or wind) is always small. As expected, the effect of discreteness is 
smaller for the outer rings. The main effect is that when discrete sources are considered, the 
very nearby sources necessary to flatten the low end of the distribution are always lacking. 
Another way of phrasing this result is that shot noise is dominated by the nearby sources. 
The absence of nearby sources had been proposed by Lezniak and Webber (1979) to explain, 
in a different context, the depletion at low grammages that was thought to be observed, 
before Webber et al. (1998) proposed a settlement to this particular issue. The total path- 
length distribution is very close to an exponential, as expected for a homogeneous source 
distribution in an infinite disk (Jones 1979). The effect of the wind is very similar to that of 
escape. It has been remarked by Jones (1978) that to a given value of (in the case where 
the wind velocity is constant in each side of the halo) can be associated an effective escape 
height L* = K/Vc = t„/2, and we have explicitly checked that the previous results apply by 
replacing L by L*. 



The previous results, as well as the results presented in our previous works, rely on 
the assumption of a constant wind presenting a discontinuity through the galactic disk. To 
probe the sensitivity of our results to this hypothesis, we consider another model in which 
the value of = Vqz varies linearly with z. This corresponds to the choice made e.g. in the 
widely used code galprop. The calculations are detailed in Appendix C, and the results 
will be used here without further justification. The path-length distribution reads 



where is the confluent hypergeometric function, also denoted iFi or M in the literature. 
The resulting path-length distribution is shown in Fig. 6 for a linear and a constant (and 
discontinuous at 2; = 0) wind. It appears that these two situations lead to very similar 
results, and one can establish a one-to-one correspondence between the parameters Vc and 
Vo of these models (see Fig. 7). It must be noted that the energy losses have not been 
considered here. Adiabatic losses are associated to the wind gradient, and their effect should 
be different for the two forms of the galactic wind: in the constant case, they are confined 
to the disk whereas in the linear case, they are present in the whole diffusive volume. 



2.4. Constant versus linear galactic wind 




with 




+ 2 . 
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3. Time-dependent diffusion equation 
3.1. Solutions 

We now turn to the problem of discreteness in time. For that, we must solve the time- 
dependent diffusion problem for an instantaneous source. Diffusion occurs independently in 
the z and r directions. Neglecting the radial boundary, pure diffusion occurs in the radial 
direction and the density can be written as 

iV(r,z,t) = ^e-^/^^*iV(.,t) 

where the function N{z,t) satisfies a time dependent diffusion equation along z. It is con- 
venient to introduce the quantities kg = K/hvariism and k^ = 2K/Vc, so that N{z,t) is a 
solution of 

dN d^N ^, . , ,dN ^, 
dm^^" ^ sign(.) — - 2kjm . 

For point-like and instantaneous sources, the radial distribution in the disk is given by (see 
Appendix A) 

N{r, z = 0,t) = exp ( ) ^ ^-^6-^'''+''-^''' sin^{A;„L} . (14) 

n=l 



AirKt V 4Kt 
where the discrete set of kn are the solutions of 

kn cota.n{knL) = —k^ — k^ (15) 

and 



L - — = L + {K + k^ 

Though it is not immediately apparent in the expression (14), the spallations are taken into 
account, through the Eq. (15) determining the Moreover, this expression is more general 
than the form that is usually found in the literature, considering only the effects of escape, 
which is obtained by replacing the right-hand side of Eq. (15) by 0. 

The result for a source that would accelerate particles for a long period of time would 
be obtained by integrating the above expression over the acceleration period. Finally, when 
energy losses are not considered, a whole energy spectrum is also simply accounted for by a 
linear superposition of delta-function sources. 
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3.2. Interpretation - Relative importance of nearby/recent sources 



The expression (14) can be written as 



N{r,z = 0,t) 



(47rA't)3/2 



1 



exp 




9{t) , 



where the function g{t) gives the correction to the purely diffusive case and takes into ac- 
count all the relevant physical effects. As such, it depends on the propagation parameters 
(spallation cross-section, galactic wind, halo height, diffusion coefficient). The importance 
of these effects depends on the position of the source in the (r, t) plane. In particular, the 
old (large t) and remote (large r) sources are more affected by all these effects. For large 
values of t, g{t) goes rapidly to zero, ensuring the convergence of the integral over time in 
Eq. (2). It also makes the function of r obtained by this integration decreases faster that 
1/r, thus ensuring the convergence of the integral over f. 

This expression will be applied in a next study to the observed distribution of sources. 
For now, we want to illustrate the possible effect of discreteness in time. For that, we 
divide the sources into several decades in age, and we compute their contribution to the 
total spectrum. The energy spectrum of a primary species can be obtained by taking into 
account the energy dependence of K in expression (14). The result is shown in Fig. 8, where 
the source distribution has been assumed to be uniform in the disk. If we denote by f the 
average distance from the sources to the Earth, the age t gives the contribution at the energy 
E for which f ~ y^I{{E)t. The more recent sources dominate the high energy tail of the 
spectrum. This is where the effect of discreteness in time is expected to be the greatest, as 
the lower decades in age contain the smallest number of sources. For the more recent decade, 
the sources have been further split into nearby (0.1 kpc < r < 1 kpc) and bulk (r > 1 kpc). 
For a rate of 3 SN explosions by century in our Galaxy, there should be about 3 nearby 
sources in the more recent age decade. It is therefore probably important to know the actual 
position and age of these sources, and to correctly model propagation from these sources, 
e.g. from Eq. 14. 



The steady-state density results from the continuous superposition of solutions for in- 
stantaneous sources, and thus can be derived from the time-dependent solution discussed 
above: 



3.3. Reformulation of the steady-state model 




— oo 



N{r, z, t) dt . 
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The integration yields 

iVstat(r, ^) = ^ exp i-h\z\) ^ c'^^Ko (r^/kl + k^^ sin(/c„L) sin {kn{L - \z\)} 

n=l 

where the Bessel function of the third kind Kq has been introduced. This expression provides 
an alternative (but is exactly equivalent) to the usual Fourier Bessel expansion over Jo 
functions. The functions Kq over which the development is performed do not oscillate, 
inducing a faster convergence. It is thus particularly well suited for sources sharply localized 
in space, as point-like sources. We have checked that this expression is fully equivalent to 
the Fourier- Bessel expansion. 

4. Summary and conclusion 

The distribution of cosmic rays sources is not continuous. The granularity of the dis- 
tribution has observable effects on the fluxes, spectra and composition, and thus should be 
considered when interpreting observed quantities. We have presented an analytical solution 
of the diffusion problem for an instantaneous point source, which takes this effect into account 
when the effects of escape through the boundaries z = ±L, convective wind and spallation 
are considered. The next step is to apply this solution to the observed local distribution of 
cosmic ray sources. 
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A. Time dependent solution of the diffusion equation 

We consider diffusion in a cylindrically symmetric box, where both disk spallations and 
galactic wind have been taken into account. The convection velocity Vc lies in the vertical 
direction and drags the particles outside so that its value is given by sign (2) Vc where a 
constant value for Vc has been assumed. The diffusion equation reads, introducing the 
quantities kg = hvariism/ K and = Vc/2K, 

1 dN (1 d f dN\ d^N] . ^ ^dN ^ 

r— + — ^ - 2fc.sign(^) — - h6iz)N . (Al) 



K dt l^r dr \ dr J dz^ J dz 

As boundaries in the radial direction r play little role in our analysis, the galactic disk 
is modelled as an infinite flat disk in the xOy plane. It is furthermore sandwiched by 
two confinement layers that extend to z = ±L. The aim of this section is to derive the 
contribution of a source 5* located at {x = 0,y = 0, z = 0} and exploding at time t = to 
the subsequent cosmic-ray density anywhere else in the Galaxy at location P {x, y, z}. The 
initial density reduces to the Dirac distribution 

Nix,y,z,t = 0) = 6{x)6{y)6{z) , (A2) 

and we would like to compute it at any time t > 0. The radial diffusion is independent 
from diffusion along z, and is not affected by any of the processes other than diffusion. As 
a result, the solution can be factorized into 

^(^'^'^'^) = 4^^^p{-4^} ^(^'^) • (A3) 



where n is a solution of 



1 Ofi O'^Ti Oti 



This preprint was prepared with the AAS I^TJt;X macros v5.2. 
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The trick is to factorize once again the time t and the vertical z behaviors so that n = 
f{z) g{t), which separates the diffusion equation into 

rv 

g' = -ag and /" - 2A;^ sign(^) /' - + -/ = (A5) 

The resulting solution may appear contrived and exceptional. Actually an infinite set of 
such functions obtains that turns out to be a natural basis for the generic solutions to 
equation (A4). The time behavior amounts to the exponential decrease g{t) = exp(— at). 
The equation on / can be solved for 2; > and z < with the appropriate boundary 
conditions f{z = ±L) = as 

f{z) = Asm\i{k{L-\z\)} when - AKa > Q 

f{z) = Asin{k{L-\z\)} when - AKa < Q (A6) 

with k'^ = a/K — /AK. The first possibility does not fulfill the disk crossing condition 
(eq. A7 with hyperbolic functions). We therefore disregard it. We then insert (A6) into 
(A5). Derivation is to be understood in the sense of distributions, because of the singularity 
of 1 2; I in z = 0. This yields 

Kf"{z) = -af{z) - 2KkA5{z) cos kL 
Inserting into (A5) gives the condition 

2/c cotan(A;L) = — fcg — k^ (A7) 

The general solution reads 

00 

N{z, t)=Y, ^-e-""* sin {kn{L - \z\)} 

n=l 

The functions sin{A;„(L — \z\)} form an orthogonal set, and it is found that 

J Sm{kn{L - \z\)}sin{kn'iL - \z\)} dz = Snn'Cn 

with 

sm2knL 

Cn = L — (A8) 

The An are found by imposing that for t = 0, the distribution is a Dirac function, 

00 

d{z) =J2'^nSm{kn{L-\z\)} 
n=l 



- 18 - 



Multiplying by sin{km{L — \z\)} and integrating over z yields 



SO that finally 

oo 

N{z,t) = ^c-^e-""*sin(A;„L) sm{K{L-\z\)} (A9) 

n=l 

and 

^^"■'"^'^^ = i^^'^^'P ("4^) sin{/c„(L- 1^1)} (AlO) 

^ ' n=l 

The radial distribution in the disk is given by 

^ ^ n=l 



B. Path- length distribution for a non homogeneous spallative disk 

This section details the derivation of expression (7), giving the grammage distribution 
in the case of an arbitrary radial distribution of spallative matter. The general method was 
sketched in Wallace (1981) to derive the cosmic ray density profile, and we present here a 
more general version which gives the grammage distribution. 

We start from Eq. (3) for the steady-state case 

1 "9 d'^Q{r,z,x)\ dg{r,z,x) g 



= K |-- irg{r,z,x)) + + g(r)5(^)5(x) - "'^^^ ' vmY^^J\r)5{z) 

with /(r) = Sism(r)/S°jjj and S^j^^ = 6.210^° cm~^ (Ferriere 1998). We perform Fourier- 
Bessel transforms using the Jq functions, 



g{r, z,x) = ^ Qiiz, x) Jo (O-^) and f{r)Q{r, z = 0, x) = ^ fi{x)Jo [Ci^ 

i=0 i=0 

with 



2 f 

Si{z, x) = ^ / pJo iCip) G{pR, z, x)dp 
/i(^) = 7277T / P-JoiCiP) f{pR)S{pR,z = 0,x) = Vaij^j(0,x) 

-'l Jo j^Q 
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where we have introduced the matrix 



2 



«ii = -j2(^ pMQp)MCjp)f{pr)dp 

These expressions are reminiscent of those of Wallace (1981) that was dedicated to a pertur- 
bative resolution of the diffusion equation in the presence of an arbitrary matter distribution 
(with no description of the path- length distribution). The generalized diffusion equation 
reads 

= K [§-,G^{^,x) + ^^1^} + qAm^) - ^vmEiJ{z) (Bl) 
The solution for z 7^ satisfying f{z = ±L, x) = is 

This expression is inserted back in the diffusion equation, taking care of the singularity of 
\z\ in which yields 

= ^.(0,x) + --tanh - ^^tanh (^^J g^^^'^^ (^2) 

The solutions of this linear set of coupled first order differential equations are 

00 

e?,(0,x) =^a,,e-^/^^e(x) (B3) 

j=0 

where the the eigenvalues of the matrix 

2KCi \ R J 

Indeed, inserting expression (B3) in (Bl), and using 

fc=0 ^ k=0 

we find an equation that can be separated in a regular part (factor of Q{x)) and a singular 
part (factor of S{x)). The former reads 

t - ^ -"'^ (f ) £ £ - « 

j=0 ^* ^ / k=0 j=0 ^ 
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In order for each coefficient of e ^^^^ to be zero, one must have 



2QKx, 



tanh I 1 2^ aikUkj = 



so that 



det 



XjOik _ , — tanh 



fc=0 







(B4) 



(B5) 



This shows that the genvalues of Aij. This equation alone is then not enough to 

compute the aij. An extra relation is provided by the singular part 



2QK 



O'jk 



j=0 k=0 



R 



which gives 



oo oo 



k=0 j=0 



ik 



(B6) 



The coefficients completely set by Eqs (B4) and (B6). 



When the galactic wind is taken into account, a more tedious derivation shows that 
Eq. (B2) should be replaced by 



= Gi{0,x) + 



2K 



\- Si coth ( 



R 



^ ' 2K 



— + Si coth ( ^ 



-1 oo 
j=0 



a 



dx 



(B7) 



so that the same results apply, provided that one makes the substitution 

V r)) '^ = H-rl' R^ 



R 



2Ci 



R 



— — tanh — - — i> \- Si coth 



and Si = 2 — + ^ 



C. Linear galactic wind in the steady-state cylindrical disk-halo model 

Resolution of the diffusion equation for a stable primary We write Vc{z) = Vqz, 
and the steady-state diffusion equation reads (see also Bloemen et al. 1993) 

K d ( ^N^ dN ( dN 1 
Developing over Bessel functions. 



K ' ' ' dz"^ K dz K ' K 
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It is convenient to rewrite this equation in a hermitic differential form, to ensure that the 
solutions form an orthogonal set of functions (see e.g. Morse & Feshbach 1953). We introduce 
rii = iVj exp(— A;^z^/2), /3 = Vq/AK and y = kz with k = ^JV^jlK, which yields, in the halo 

n/' - rii [ai + y^) = 

where 

The solutions are of the form, taking into account the condition Uii^z = ±L) = 0, 




ni = Bie " ' {(t){ — -; k^z^ ] - - )^ , f ( — —, -; 



4 '2' / L , /3 + ai 3 ,2.2 \ V 4 '2 



4 '2' 



where (/> is the confluent hypergeometric function, also noted iFi. The value of Bi is found 
by integrating the diffusion equation through the disk, so that 



2 dui 
dz 



This gives Bi = 2hqi/Ai with 



^^ = ^4t^:-^ f + -^sL (Ci) 



3^ 3 2^2 

4 ' 2' 



The final solution is thus obtained as 




4 '2'^^ 



^ ' ^ ^ A: V^*i?y r V 4 '2' / L , /3 + ai 3 ^ V 4 '2' 

(C2) 



The density in the disk is thus given by 

2hqi 



iV(r,. = 0) = 5:^^Jo(c^) (C3) 

i 

It can be shown that this expression reduces to the usual expressions in the case of a vanishing 
wind. 
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The path-length distribution The dependence in a is very simple and the path-length 
distribution is obtained by inverse Laplace transform as 



N{r, z,x) = Y, '^hqiJo (Ci^) exp (^-^^ e(x) 



with 



5 K Q 3 V^L 



_ mvEL^ V4 2K)i?2' 2' 2K 
~ 2K ^ /3 K Cf 1 VqL' 



4 2K)i?2' 2' 2/s: 



D. Remark about the time evolution of the mean grammage 

The mean grammage of the CR emitted by a single source can be expressed from the 
distribution Q (r , x, t) as 



xit) = ^y.. (Dl) 



IJ 




fOO 

dxxQ{f, X, t) 

) 


L 




poo 

/ dxQ{f,x,t) 
'o 



In this expression, the averaging process is understood to be performed over the whole spatial 
distribution of cosmic rays. The time derivative of this expression yields 



dX 




/oo r r r 

dx X dQ{f,x,t)/dt / / / d^f / dx dQ{f,x,t)/dt 

- X -^-^i,, ,oo (D2) 




/•oo 

d^r I dxQ{f, X, t) 
Jo 




d f / dxQ{f, X, t) 



The denominator is simply N{t), the total number of CR in the diffusive volume. Using the 
fact that 

— — = —mpcmr) 
dt w 

where n(r) is the interstellar gas density and integrating the first term by parts, we finally 
find, 

dX ^ , , s, 1 dN 

— = m^crisnit) - Xit) -— (D3) 

where nsn{t) is defined as in (8). The last term is missing in Eq. 4.3 of Higdon and Lingenfelter 
(2003). This term is positive and represents the change in grammage due to the escape of a 
fraction of cosmic rays between times t and t + dt. 
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As stressed in the text, a quantity which has a greater physical importance to us is the 
mean grammage of the CR that reach the Earth at time t, i.e. 



dxxQ{f, X, t) 

x{r, t) = ^ (D4) 

dxQ{f, X, t) 



Following the same procedure as above, its time evolution is given by 

dxiv t\ d 

— — = fh(3cn{f^ — x{r,t) — {\Yiw{r,t)} (D5) 
where the definition of w{r^t) is given in (10). 
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qrammaqe distribution dN/dx, loq scale 
-10.2 




Fig. 1. — The upper curves represent the path- length distribution (x in g cm~^) for L = 4 
kpc and a homogeneous disk distribution corresponding to the local matter density, and with 
a realistic gas distribution, including the H2 ring, radially distributed according to Ferriere 
(1998). The lower curves are obtained by taking into account the spallation cross-section 
of Oxygen. The feature at x ~ 2 g cm~^, visible in the Higdon points (dotted line) is not 
reproduced by the analytical result. 
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Fig. 2. — Secondary-to-primary ratio B/C, as a function of the diffusion coefficient normali- 
sation Kq in kpc^/Myr, for a homogeneous matter distribution (sohd hne) and for a reahstic 
gas distribution (dashed line). 
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Fig. 3. — Path-length distributions (x in g cm~^) are shown for a point source located at a 
distance r = 1 kpc, and for different propagation conditions. The thick line corresponds to 
free diffusion (no boundary, no wind). The thin dotted lines show the influence of escape 
for L = 1, 2 and 5 kpc. The higher grammages, corresponding to longer paths, are more 
suppressed by escape, which shifts the maximum of the path-length distribution towards low 
grammages. The effect of the convective wind, displayed as solid lines for = 1 and 2 kpc, 
is similar though quantitatively different for low grammages (see text for a discussion). The 
overall normalisation is decreased when escape and/or wind are considered, which reflects 
the corresponding decrease of flux. 
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Fig. 4. — The path-length distribution [x in g cm~^) is shown for a point source and for 
different r. The two situations L = 5 kpc and L infinite are shown for comparison. The 
short grammages are depleted for remote sources, as the probability to travel a long distance 
without crossing the disc is small. This effect is all more pronounced that L is small. The 
normalisation also decreases with L, which reflects the corresponding decrease of flux. 
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grammage distribution dN/dx, log scale 




2 4 6 8 '0 12 ^4 16 18 20 



X (g/ cm^) 

Fig. 5. — The contributions to the path-length distribution (thick hne) from continuous 
sources distributed in several rings (thin lines) are displayed for L = 5 kpc. The path-length 
distribution from a set of 100 point sources drawn randomly in the central ring (inner 1 kpc) 
is also shown (dotted line). The effect of discreteness is most visible for low grammages. For 
the other rings, the difference between the continuous distribution and 200 discrete sources 
is of the order of a percent. The effect of galactic wind is very similar to that of L, replacing 
LhyL* = K/V, = 
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Fig. 6. — Path-length distributions for a hnear (symbols, Vq = 0.018 kpc Myr~^ kpc~^) and 
a constant galactic wind (line, = 0.01 kpc Myr~^). It appears that the effect of a linear 
or constant wind are very similar. For each value of Vq, it is possible to find a value of Vc 
yielding a grammage distribution which is very similar, with an accuracy better that one 
percent. 
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Fig. 7. — Correspondence between the values of Vc for a constant wind {x axis) and Vq for a 
linear wind {y axis) giving approximately the same path-length distributions, for two values 
of L. 
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Spectrum dN/dR 
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Fig. 8. — The contribution to the spectrum from different age classes is displayed for a 
uniform source distribution. In the more recent age class (0.001 Myr < t < 0.01 Myr), the 
contribution from nearby sources (r < 1 kpc) has been separated from the rest (r > 1 kpc). 
The effect of time discreteness will be more important in this recent age class, which contains 
less sources. 



